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CHAPTE"  1 


INTRODUCTION 

This  present  work  is  part  of  a  study  of  the  behavior  of 
high-strength  concrete  under  monotonic  biaxial  compressive  loading. 

A  model  of  concrete,  consisting  of  a  square  mortar  plate  with  nine 
coarse  aggregate  circular  inclusions  (see  Fig.  2.1),  is  analyzed 
using  the  Finite  Element  Method.  The  analytical  results  are  then 
compared  with  results  of  experimental  tests  of  the  same  model. 

The  analysis  takes  into  account  the  nonlinear  behavior  of  the 
mortar  using  the  constitutive  equations  proposed  in  Ref.  1.  These 
constitutive  equations  have  proved  adequate  for  the  prediction  of 
the  behavior  of  concrete  under  biaxial  states  of  stress. 

The  significance  of  the  bond  between  the  coarse  aggregate  and 
the  mortar  is  also  studied  using  an  interface  element  developed  in 
this  work.  Relative  displacement  of  the  two  materials  occurs  if  the 
strength  of  the  bond  is  exceeded. 

Chapter  2  gives  a  brief  account  of  the  behavior  of  high- 
strength  concrete  and  presents  a  comparison  with  normal -strength  con¬ 
crete.  The  properties  of  the  materials  used  in  the  analysis  are  also 
presented.  The  constitutive  equations  used  for  the  mortar  are 
described  in  detail  in  Chapter  3  and  the  linear  and  nonlinear  finite 
element  analysis  of  the  concrete  model  are  developed  in  Chapter  4. 
Chapter  5  shows  the  modeling  of  the  interface  between  the  mortar  and 
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the  coarse  aggregate  and  the  bond  properties  between  the  two  mater¬ 
ials.  Chapter  6  gives  the  analytical  results  and  compares  them  with 
experimental  ones.  Finally,  some  conclusions  about  the  important 
factors  affecting  the  stiffness  and  strength  of  the  concrete  model 
are  listed  in  Chapter  7. 


CHAPTER  2 


HIGH-STRENGTH  CONCRETE 


2. 1  Introduction 

The  concrete  model  shown  in  Fig.  2.1  was  first  proposed  in 
Ref.  2  for  the  study  of  the  behavior  of  normal-strength  concrete 
subjected  to  biaxial  loading.  In  the  present  work,  the  same  model 
is  used  to  study  the  behavior  of  high-strength  concrete  under  biaxial 
loads.  Thus,  this  chapter  gives  a  brief  presentation  of  applications 
and  properties  of  high-strength  concrete  and  summarizes  the  most 
important  differences  in  behavior  between  high-strength  concrete  and 
normal-strength  concrete.  The  properties  of  the  component  materials 
in  the  concrete  model  are  also  presented. 

2.2  Applications  of  High-Strength  Concrete 

In  recent  years,  efforts  have  been  made  to  improve  the  com¬ 
pressive  strength  of  concrete.  Nowadays,  high-strength  concrete  is 
becoming  increasingly  common.  Among  other  applications,  it  has  been 
used  in  high-rise  buildings  where  oversized  columns  can  be  avoided 
in  the  lower  floors,  or  in  prestressed  flexural  members  permitting 
larger  values  of  span-to-depth  ratio.  It  also  has  great  potential 
of  use  in  structures  in  which  the  concrete  is  subjected  to  a  biaxial 
state  of  stress  such  as  large  shells,  containment  vessels  and  tunnels. 
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MORTAR 

_ J _ 


AGGREGATE 


Fig.  2.1  The  concrete  model 
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2.3  Properties  of  High-Strength  Concrete 

2.3.1  Uniaxial  Compressive  Strength.  Although  there  is  no wei  1 
defined  boundary  between  normal -strength  and  high-strength  con¬ 
crete,  some  authors  (see  Ref.  3)  classify  arbitrarily  concrete  with 

a  specified  cylinder  compressive  strength  of  at  least  6,000  psi  as 
being  high-strength  concrete.  Concretes  with  specified  cylinder 
compressive  strength  in  the  range  of  8,000  to  10,000  psi  have  been 
used  successfully  with  conventional  technology  and  materials  but  with 
careful  selection,  proportioning  and  quality  control.  In  this  study, 
the  uniaxial  compressive  strength  of  the  concrete  model  ranged  between 
6,000  and  8,000  psi. 

2.3.2  Stress-Strain  Relations.  Plain  concrete  has  some 
amount  of  ductility.  This  ductility,  however,  decreases  with  increas¬ 
ing  concrete  strength.  The  stress-strain  relation  up  to  ultimate 
strength  becomes  almost  a  straight  line  as  the  concrete  strength 
increases  (see  Fig.  2.2  and  Ref.  3).  Note  that  there  is  a  descending 
branch  in  each  curve  after  the  maximum  stress  has  been  reached.  Also, 
the  maximum  strain  at  failure  in  compression  is  lower  at  higher  con¬ 
crete  strengths.  The  maximum  ultimate  strain  may  be  below  0.003  for 
higher-strength  concretes.  As  it  may  be  seen  in  Fig.  2.2,  the  modulus 
of  elasticity  is  greater  for  higher  strength  concrete. 

2.3.3  Microcrackinq.  The  differences  in  behavior  between 
high-strength  and  normal-strength  concrete  as  shown  in  Fig.  2. 2  may 
be  explained  by  differences  in  microcracking.  For  higher-strength 
concretes  there  is  less  cracking  at  the  interface  between  the  aggregate 
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Fig.  2.2  Stress-strain  curves  for  concrete  (Ref.  3) 

and  the  mortar  than  for  lower- strength  concretes  and  the  stress- 
strength  ratio  at  which  microcracks  begin  to  form  continuous  crack 
patterns  is  higher.  Therefore,  the  stress-strain  curve  is  steeper 
and  closely  linear  up  to  a  higher  stress-strength  ratio  (see  Ref.  4), 
and  the  number  of  continuous  crack  paths  is  smaller  for  higher-strength 
concrete  resulting  in  a  decrease  in  the  redundancy  present  in  the 
material.  This  is  an  explanation  for  its  lack  of  ductility. 

2.3.4  Failure  Mode.  High-strength  concrete  behaves  more 
like  a  homogenous  material  than  normal-strength  concrete  and  therefore 
their  failure  modes  in  uniaxial  compression  are  different.  The  frac¬ 
ture  surface  of  normal-strength  concrete  generally  follows  the  con¬ 
tour  of  the  coarse  aggregate  surface  in  inclined  planes  forming  a 
cone  of  rupture.  In  high-strength  concrete,  failure  occurs  in  a  plane 
parallel  to  the  applied  load  and  passing  through  the  aggregate  and 
the  mortar  (see  Ref.  4). 
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2.4  Component  Materials  Used  in  the  Concrete  Model 

Typical  stress-strain  curves  for  the  mortar  and  the  coarse 
aggregate  used  in  the  concrete  model  are  shown  in  Fig.  2.3.  As  it 
may  be  seen,  the  coarse  aggregate  behaves  almost  linearly  up  to 
ultimate  strength  and  the  mortar  shows  some  nonlinearities  for  loads 
higher  than  approximately  40  percent  of  ultimate  strength.  Thus, 
in  this  study,  the  coarse  aggregate  will  be  assumed  to  be  elastic  and 
the  nonlinear  behavior  of  mortar  will  be  taken  into  account. 

The  coarse  aggregate  and  the  mortar  used  in  this  study  were 
obtained,  respectively,  from  limestone  rock  and  from  a  mix  of  natural 
sand  and  type  I  cement  (w/c  =  0.35,  s/c  =  2.0).  The  elastic  proper¬ 
ties  of  the  two  materials  obtained  from  the  average  of  three  tests 
on  cylinder  specimens  are  presented  in  Table  2.1. 


TABLE  2.1  Elastic  Properties  of  the  Mortar  and  the  Coarse  Aggregate 

Material  Coarse 

Property  Aggregate  Mortar 

Compressive 

Strength  (psi)  12,850  10,970 

Tensile 

Strength  (psi)  —  1,029 

Modulus  of  c  c 

Elasticity  (psi)  4.49  x  10°  6.74  x  10° 


Poisson's 

Ratio 


0.058 


0.25 


CHAPTER  3 


DESCRIPTION  OF  A  CONSTITUTIVE  MODEL 


3.1  Introduction 

The  behavior  of  the  coarse  aggregate  and  the  mortar  under 
short-term  monothonic  compression  was  discussed  in  detail  in  Chapter  2. 
In  uniaxial  compression  the  stress-strain  curve  for  the  mortar  may  be 
assumed  linear  for  levels  of  load  up  to  40-50  percent  of  ultimate 
strength.  Beyond  this  point,  inelastic  behavior  is  clearly  obtained 
and  must  be  considered  in  the  analysis. 

In  this  chapter  a  model  is  presented  applicable  to  the  descrip¬ 
tion  of  the  behavior  of  the  mortar  under  multiaxial  states  of  stress. 
The  model  parameters  are  estimated  from  experimental  results.  As 
outlined  below,  the  development  of  the  model  follows  the  formalism 
of  the  classical  theory  of  plasticity.  The  model  was  first  proposed 
in  Ref.  1  in  which  details  may  be  found. 

3.2  Description  of  the  Model 

Following  the  theory  of  plasticity,  the  increment  of  strain 
dc - .  is  taken  to  be  the  sum  of  an  increment  of  strain  resulting  from 

'  J 

elastic  behavior  de?j  and  an  increment  of  strain  resulting  from  inel¬ 
astic  (plastic) behavior  dc?.. 

*  sJ 


* 


(3.1) 
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Assuming  that  the  elasticity  of  the  material  is  isotropic,  the 

e 

elastic  increment  of  strain  d  is  given  by: 


e  e 

d  eij  =  C i j k £  dak£ 


(3.2) 


where  C i j k £ 
by: 


is  the  isotropic  elastic  compliance  tensor  and  is  given 


Cijk£  2G  6ik  5j£  +  '  6^  ^  6ij  6  k£^  (3’3) 

(Indicial  notation  is  conveniently  used.  Indices  assume  the  values 
1,  2,  3.  The  Kronecker  delta  6*.-  is  defined  as  6 .  ■  =  0,  if  i  f  j 

'  J  *  J 

and  <$..  =  1,  if  i  *  j).  K  and  G  are  the  elastic  bulk  and  shear 

*  J 

moduli  defined  as: 


1(1-  2  u) 


(3.4) 


G 


E 

2  (1  +  u) 


(3.5) 


E  and  u  are  the  Young's  modulus  and  the  Poisson's  ratio  respectively. 
An  alternative  form  of  Eq.  (3.2)  is: 


(3.6a) 
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d  e6 

_  da 
'  K 

(3.6b) 

=  dee. 
ij 

+  5  “ 

(3.7) 

“  dsij 

*  do5ij 

(3.8) 

0  p  e 

d  e  is  the  elastic  increment  of  deviatoric  strain,  d  e  =  d  e  ^ 
is  the  elastic  increment  of  volumetric  strain,  do  =  d  o„  is  the 
increment  of  hydrostatic  stress  and  d  S-.  is  the  increment  of  devia- 

'  J 

toric  stress. 

An  important  assumption  in  the  calculation  of  the  plastic  incre¬ 
ment  of  strain  d  £?-,  following  the  theory  of  plasticity,  is  the  exis- 

'  vi 

tence  of  a  yield  function.  In  the  present  model,  the  arguments  of  the 
yield  function,  F,  are  the  state  of  stress,  a^.,  and  a  parameter,  k, 
that  reflects  the  history  of  plastic  deformation: 

F  (  a.y k)  =  0.  (3.9) 

In  the  space  of  stresses  o. the  yield  function  may  be  repre- 

*  J 

sented  by  a  surface.  For  plastic  deformation  to  occur  the  material  must 
be  subjected  to  a  state  of  stress  that  lies  on  this  surface,  F  =  0 
(yield  surface).  If  the  state  of  stress  is  in  the  interior  of  the 
convex  region  bounded  by  the  surface,  F  <  0,  only  elastic  deformation 
occurs.  The  value  of  the  parameter  k  changes  so  that  the  state  of 
stress  satisfies  F  =  0  during  plastic  deformation.  Thus,  the  region 


1 


L 
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defined  by  F  >  0  is  the  set  of  states  of  stress  which  cannot  be  ob¬ 
tained  without  further  plastic  deformation.  In  this  way,  an  increment 
of  stress  do^j  can  be  considered  as  loadino,  if  it  is  directed  towards 
the  exterior  of  the  convex  region  bounded  by  the  yield  surface,  as 
unloading  if  it  is  directed  towards  the  interior  of  the  convex  region 
and  as  neutral  loading,  if  it  is  tangent  to  the  yield  surface: 

1)  Loading: 

F  (o4.»k)  =  0  and  do..  >  0  (3.10a) 

1J  30  ij  1J 

2)  Unloading: 

F  (o.,,k)  =  0  and  -3-^-  do..  <  0  (3.10b) 

Ij  dC-j  1J 

3)  Neutral  Loading: 

F  (0ij,k)  =  0  and  do..  .  0  (3.10c) 

(It  is  assumed  that  F  is  continuously  differentiable.)  Plastic  deform¬ 
ation  occurs  only  during  loading.  The  yield  function  F  is  assumed  iso¬ 
tropic,  involving  invariants  of  the  stress  tensor  o^,  namely,  the 
hydrostatic  stress  (mean  normal  stress) 


and  the  shearing  stress  intensity 


13 


/!/6  '[(oj  -  o^)2  +  (cr2  -  o3)2  +  (03  -  Qj)2] 


1/2 


(S..  =  o-  -  -  a  5...  a,,  a?,  a,  are  the  principal  values  of  the  stress 

1 J  I J  •  J  1  £  J 

tensor.)  It  is  taken  of  the  general  form 


F  (a,  T,  k)  =  f  (a,  T)  -  k 


(3.11) 


Upon  loading,  F  (a^j.k)  =  0  and  F  (a^  +  d  a^.,  k  +  d  k)  =  0. 
Thus, 


dF-;577  d* 


d  0  .  .  -  d  k  =  0 

30 ij 


therefore 


d  k  =  — •  d  a-- 

30^-  ij 


The  plastic  increment  of  strain  d  £.?  is  written  as 

*  \J 


d  P  P  =  de'3  +  —  dc^ii 
o  Cij  a  e-j  +  3  e  0^ 


(3.12) 


(3.13) 


with  d  e^j  (the  plastic  increment  of  deviatoric  strain)  and  d  e*3  (the 
plastic  increment  of  volumetric  strain)  given  by: 


d  e  P  =  3.  d  k 

a  eij  H  a  K  2  T 


(3.14a) 


I 
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d  eP  =  l  d  k  3 


(3.14b) 


with  H  and  3,  in  general,  functions  of  o  and  T  . 

In  order  to  describe  the  physical  meaning  of  H  and  3,  it  is 
convenient  to  define  the  intensity  of  the  increment  of  inelastic 
shear  distortion  d  ep  as: 


d'  -  t2  d  ejt  d  ef.  ] 


1/2 


"k«. 


(3.15) 


Using  Eq.  (3.14a),  it  may  be  seen  that,  during  loading,  H 
must  be  positive,  since  d  k  >  0,  T  >  0  and  d  e^  is  taken  in  the  same 
direction  as  ..  It  is  easily  obtained  that 


d  e 


P 

H 


2 

d  k 
H 


\a 


Since  T 


2 


1 

2 


S ^  Sk  ,  it  follows  that 


2  d  e 


P 

k£ 


d  k 
H 


2 


and 


c  2  d  eu  d  ekP,  ] 


1/2 


d  k 
H  • 


Final  ly 


Since  d  k  may  be  understood  as  an  increment  of  loading  and  dep  is  the 
intensity  of  the  increment  of  plastic  shear  distortion,  H  may  be 
interpreted  as  a  plastic  shear  modulus. 

Combining  equations  (3.14b)  and  (3.16),  it  follows  that 

d  eP  =  Bd  ep.  (3.17) 

Thus  B  may  be  understood  as  the  ratio  of  inelastic  volumetric  deforma¬ 
tion  to  inelastic  shear  deformation.  It  is  referred  to  as  the  inelas¬ 
tic  dilatancy  factor.  Since  d  ep  >  0,  B<0  means  inelastic  contraction 
and  3>0  means  inelastic  dilatancy. 

The  yield  function  F  must  exhibit  what  is  sometimes  referred 
to  as  pressure  sensitivity  of  inelastic  behavior.  Thus,  it  is  assumed 
that  as  magnitude  of  the  hydrostatic  stress  (a|  increases,  the  shearing 
stress  intensity  for  which  inelastic  behavior  may  occur,  also  increases. 
A  typical  yield  surface  in  a-T  space  is  shown  in  Fig.  3.1.  The  para¬ 
meter  k  is  also  known  as  the  hardening  parameter  and,  if  the  state  of 
stress  is  on  the  yield  surface,  it  is  related  to  o  and  T  through  the 
equation  F  =  0  (Eq.  (3.11)). 

Thus,  H  and  B  may  be  taken  as  functions  of  (o,T),  [a,  k)  or 
(t,  k)  equivalently.  However,  since  ultimate  strength  of  the  material 
is  obtained  when  k  reaches  a  limit  value,  it  is  more  convenient  to  take 


Fig.  3.1.  Typical  yield  surface  in  a  -  T  space 


H  and  8  as  functions  of  (T,k)  or  (a,k)  rather  than  (o,T)  since  k 
provides  a  measure  of  the  proximity  of  the  state  of  the  material  to 
the  ultimate  strength. 

3.3  Constitutive  Equations 

Eqs.  (3.14)  may  be  rewritten  as: 
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The  plastic  compliance  tensor  C  is  obtained  using  Eq.  (3.12): 


d  k  = 


3akt 

dak£j 

3  f 

3SW 

dSk£ 

3  f  . 

+  t —  d  a 
3o 

3  f 

hi  d 

3  f 

C  + 

Ft 

2  T  d 

k£  3  a 

Since 


Sksi  CT  U  "  o6ka 


d  Ske  =  dok£  ’  do  6k£ 


it  is  seen  that 


3  f  Sk£ 

d  k  Ft  2~T  d(Jkil 


3  f  Sk£  do  5.  + 

5  T  2  T  ^  9  a 


However, 


3  f  sk£ 

ST  2T  dc\i 
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d  a 


3  d  V  6kZ  * 


therefore, 


d  k  = 


H  llii  +  I  H  6kz 

3T  It  3  3  o 


do 


kr 


It  follows  from  Eqs.  (3.13)  and  (3.14)  that: 


C  p  -  - 
i  jk£  H 


Sii  1 
2~T  +  3  6  6ij 


!!  , 

9T  2T  3  9o  k£ 


(3.19) 


In  general.  cj^j. 


Finally,  the  incremental  constitutive  equations  corresponding 
to  loading  may  be  written  as 


dcij  =  Cijk£  dGk£ 


(3.20) 


with 


f  =  C  e  +  r  P 
uijk£  Lijk£  S'jk£  * 


Again,  in  general,  Cijka  f  C^.. 

In  the  cases  of  unloading  and  neutral  loading,  the  compliance  tensor 
Cijk£  is  e9ua1  t0  the  e^astic  compliance  tensor  C..®^  and  the  increment 
of  strain  d  e^  is  given  by: 
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d  eij  Ci jkn  dok2,  *  (3'21^ 

3.4  Parameter  Estimation 

The  yield  function  F  has  been  taken  of  the  general  form: 

F  (a,  T,  k)  =  f  (a,  T)  -  k  (3.11) 

A  simple  special  case  of  Eq.  (3.11)  is: 

F  (a ,  T,  k)  =  Al0  +  A2T  +  A3  -  k  =0 

Using  the  following  conditions, 

1)  a ^  =  - 1  f c  I  «  °2  =  °3  *  13  dnd  **  =  1 

2)  =  ©2  =  -1.16  | f 1 c | ;  =  0  and  k  =  1 

3)  Oj  =  -0.40  | f ^ |  S  °2  =  0  3  =  0  and  k  =0, 

the  values  of  A^,  A2  and  A^  may  be  found.  Thus 

F  (0  ,  T,  k)  3  0.69  -pp-j-  +  3.285  ypy  -.0667  -  k 


(f^  is  the  strength  of  the  material  in  uniaxial  compression.) 
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Conditions  1,  2  and  3  represent* respectively* 

1)  ultimate  strength  in  uniaxial  compression. 

2)  ultimate  strength  in  equal  biaxial  compression.  (The 
factor  1.16  allows  for  the  increase  in  strength  under  combined  com¬ 
pressive  stresses. ) 

3)  Initial  yielding  in  uniaxial  compression. 

The  parameter  k  can  assume  values  between  0  (initial  yielding) 
and  1  (ultimate  strength).  Thus,  it  can  be  used  as  an  indicator  of 
the  proximity  of  the  yield  surface  to  the  material  ultimate  strength. 

A  graphical  representation  of  the  yield  function  (3.22)  in  -  o ^ 
space  and  in  a  -  T  space  is  shown  in  Figs.  3.2  and  3.3  respectively. 

The  plastic  shear  modulus  H  will  be  taken  to  be  a  function  of 
the  magnitude  of  the  hydrostatic  stress  |o|  and  the  hardening  parameter 
k.  As  shown  in  Figs.  3.1  and  3.3,  the  yield  function  used  in  this 
study  is  such  that  the  shearing  stress  intensity  T  increases  as  the 
magnitude  of  the  hydrostatic  stress  |o|  increases.  Since  more  inelas¬ 
tic  shear  deformation  is  obtained  under  higher  shear  loading  (higher 
T),  is  follows  that  H  must  be  a  decreasing  function  of  T  and  hence  of 

M* 

Using  Eq.  (3.16)  and  calculating  the  values  of  d  k  and  dep  from 
stress-strain  curves  measured  in  uniaxial  compression  tests,  the 
values  of  H  shown  in  Table  3.1  were  obtained.  Within  each  interval, 

H  is  assumed  to  vary  as  follows: 

1 

a  +  b  k 


H 


Table  3.1  Relationship 

between  k  and  H  in  uniaxial  compression 

k 

H  (psi) 

0 

16,000 

0.35 

4,870 

0.50 

2,800 

0.70 

1,600 

1.00 

575 

The  values  of  a  and  b  may  be  calculated  for  each  interval  yielding 
the  expressions  in  Table  3.2. 


Table  3.2  Expressions  for  H  as  a  function  of  k  in 
uniaxial  compression. 


_ k _ 

0  ^  k  <  0.35 
0.35  5  k  5  0.50 
0.50  i  k  <  0.70 
0.70  5  k  5  1.00 


_ H  (psi) _ 

H  =  1/(4.081  x  10'4k  +  6.250  x  10"5) 
H  =  1/(1.012  x  10~3k  -  1.489  x  10‘4) 
H  =  1/(1.339  x  10'3k  -  3.125  x  10“4) 
H  =  1/(3.714  x  10'3k  -  1.975  x  10"3) 


In  uniaxial  compression  (  <  0,  02  =  aj  ~  0); 


o 


l 

3 


al 


and 


I 
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1 

rr 


Using  the  yield  function  given  by  Eq.  (3.22),  it  is  seen  that,  beyond 
the  initial  yielding. 


2  +  3  k 


ifcl 


and  therefore 


2  +  3  k 
15 


or 


2  +  3  k 

~  I  a  F 


In  order  to  generalize  the  results  given  in  Table  3.2  and  obtain 
expressions  for  H  valid  for  any  state  of  stress,  the  values  of  H  in 
uniaxial  compression  were  divided  by  15  and  multiplied  by 


2  +  3  k 

I.  Q__| 


Thus,  the  expressions  for  H  shown  in  Table  3.3  were  determined. 


24 


Table  3.3  Expressions  for  H  under  multiaxial  states 
of  stress  (C  =  {2  +  3  k)/]o  /f^l). 


_ k _ 

0  *  k  -  0.35 
0.35  *  k  *  0.50 
0.50  *  k  *  0.70 
0.70  *  k  *  1.0 


_ H..(psi) _ 

H  =  C/(6. 122  x  10'3k  +  9.375  x  10"4) 
H  =  C/ (1.518  x  10'2k  -  2.234  x  10'3) 
H  =  C/(2.009  x  10“2k  -  4.688  x  10'3) 
H  *  C/ ( 5 . 571  x  10"2k  -  2.963  x  10'2) 


The  inelastic  dilatancy  factor  8  will  be  taken  to  be  a  function 
of  the  parameter  k  only.  A  more  elaborate  analysis  could  include  the 
effect  of  the  hydrostatic  stress  a  in  the  expressions  for  6. 

Using  £q.  (3.17)  and  calculating  the  value  of  dep  and  dep  from 
stress-strain  curves  measured  in  uniaxial  compression  tests,  the  value 
of  6  shown  in  Table  3.4  were  obtained. 

In  order  to  check  the  accuracy  of  the  expressions  for  H  and  B 
shown  in  Tables  3.2  and  3.4  respectively,  the  stress-strain  curve  in 
uniaxial  compression  used  to  obtain  the  two  parameters  was  reproduced 
and  the  comparison  is  shown  in  Fig.  3.4. 


1 
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Table  3.4  Relationship  between  k  and  8  in  uniaxial 
compression. 

k  8 

0  ^  k  *  0.135  6  =  -49.60k  +  2.92 

0.135  <  k  *  0.219  8  =  66.31k  -  12.73 

0.219  0.354  8  =  -14.44k  +  4.95 

0.354  S  k  S  0.521  8  =  -  3.05k  +  0.92 

0.521  S  k  £  0.615  8  =  8.09k  -  4.88 

0.615  *  k  S  0.708  8  =  -9.89k  +  6.17 

0.708  -<k^  0.792  8  =  15.12k  -  11.53 

0.792  S  k  S  0.875  8  =  -11.81k  +  9.79 

0.875  *  k  S  i.o  8  =  8.32k  -  7.82 


CHAPTER  4 


THE  FINITE  ELEMENT  ANALYSIS  OF  THE  CONCRETE  MODEL 

4.1  Introduction 

The  concrete  model  to  be  analyzed  in  this  study  is  shown  in 
Figure  2.1.  The  coarse  aggregate  will  be  assumed  to  be  elastic  up 
to  ultimate  strength  since  test  results  in  uniaxial  compression  showed 
that  this  is  true  to  a  very  good  approximation.  Inelastic  deformation 
in  the  mortar  will  be  taken  into  account  using  the  constitutive  model 
described  in  Chapter  3. 

As  the  model  and  the  loading  are  symmetric  about  the  two 
orthogonal  axes,  it  is  necessary  to  analyze  one  quadrant  only.  A 
discretization  of  this  quadrant  by  finite  elements  is  shown  in  Fig. 
4.1.  As  it  may  be  seen,  a  six-node  element  is  used  to  represent  the 
coarse  aggregate  and  an  eight-node  element  is  used  to  represent  the 
mortar. 

The  objective  of  this  chapter  is  to  describe  the  linear  and 
nonlinear  finite  element  analysis  of  the  concrete  model. 

4.2  Linear  Analysis 

This  section  describes  the  analysis  of  the  concrete  model 
assuming  both  materials  (aggregate  and  mortar)  to  be  isotropic  and 
elastic.  Section  4.3  will  describe  the  nonlinear  analysis,  i.e.,  the 
analysis  accounting  for  the  nonlinear  behavior  of  the  mortar. 
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The  linear  analysis  may  be  described  briefly  as  follows: 

Step  1.  Given  the  elastic  properties  of  the  materials,  a 
relationship  between  nodal  forces  and  nodal  displacements  for  each 
element  is  found: 

F1  =  K1  U1 

/\y 

where 


[xj. 


X1,  Y1'  ]T 

m 5  m  J 


[ o i ,  v^ ,  ^2’  v2» 


m’ 


The  superscript  T  indicates  that  the  transpose  of  the  superscripted 
matrix  must  be  taken.  X1,  and  y’.  are  the  forces,  in  the  x  and  y  dir- 

J  J 

ections,  respectively,  applied  at  node  j  of  element  i.  u.  and  v. 

J 

are  the  displacements,  in  the  x  and  y  directions,  respectively,  of  node 
j.  m  is  the  number  of  nodes  of  the  element.  K1  is  the  stiffness 
matrix  of  element  i. 

Step  2.  The  element  stiffness  matrix  K1  and  the  load  vector 

A-/ 

F1  are  modified  for  the  support  conditions. 

Step  3.  A  relationship  between  nodal  forces  and  nodal  dis¬ 
placements  for  the  structure  is  found,  i.e.,  the  structure  stiffness 
matrix  is  assembled  from  the  element  stiffness  matrices: 

=  K  U 

A/  A/ 


\ 


F 
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L  =  ^Xl’  Yl»  X2’  Y2»  •  •  •  ‘  XP»  Yp^ 

^  =  ^2*  v2*  •  •  •  » ^p* 

X.  and  Y.  are  the  forces,  in  x  and  y  directions  respectively,  applied 
J  J 

at  node  j,  pis  the  number  of  nodes  in  the  structure.  K  is  the 

A/ 

stiffness  matrix  of  the  structure. 

Step  4.  The  system  of  equations  for  the  displacements  is  solved 
using  Gauss  elimination.  Formally: 

U  =  [K]'1  F 

Step  5.  After  the  displacements  are  determined,  it  is  possible 
to  find  strains  and,  therefore,  stresses  at  any  point  in  the  structure. 

In  each  element,  some  points  also  used  in  the  calculation  of  the  element 
stiffness  matrix  (integration  points)  are  chosen. 

4.2.1  Derivation  of  the  Element  Stiffness  Matrix.  The  finite 
elements  used  in  this  study  are  called  isoparametric.  The  same  inter¬ 
polation  functions  used  to  relate  the  coordinates  of  any  point,  within 
the  element,  to  the  coordinates  of  the  element  nodes  are  also  used  to 
relate  the  displacements  of  any  point,  within  the  element,  to  the 
displacements  of  the  element  nodes.  Thus: 

m 

x  =  l  N  •  x .  (4.1a) 

j=l  3  3 


L 
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m 

y  =  Z  N.  y.  (4.1b) 

j=l  J  3 

m 

u  =  Z  N.  u-  (4.2a) 

j=l  J  3 

m 

v  -  £  N.  >/.  (4.2b) 

j=l  3  3 

x,  y  are  the  coordinates  of  any  point  within  the  element,  x^,  y.  are 

J  J 

the  coordinates  of  node  j.  u,  v,  are  the  displacements  of  any  point 
within  the  element,  u-,  v-  are  the  displacements  of  node  j.  N,  is 

J  J  J 

the  interpolation  function  corresponding  to  node  j. 

An  interpolation  function  assumes  a  value  equal  to  one  when 
evaluated  at  its  corresponding  node  and  a  value  equal  to  zero  when 
evaluated  at  any  other  node.  For  example,  the  interpolation  functions 
used  for  the  eight-node  element  (see  Fig.  4.2)  are  given  by 

Nx  «-l/4(r  +  s  +  1)  (1  -  s)  (1  -  r) 

N2  =  l/4(r  -  s  -  1)  (1  -  s)  (1  +  r) 

N3  =  l/4(r  +  s  -  1)  (1  +  s)  (1  +  r) 

N4  =-l/4(r  -  s  +  1)  (1  +  s)  (1  -  r) 


N5  =  1/2(1  -  s)  (1  +  r)  (1  -  r) 
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N6  =  1/2(1  -  s)  (1  +  s)  (1  +  r) 

N7  =  1/2(1  +  s)  (1  +  r)  (1  -  r) 

Ng  =  1/2(1  +  s)  (1  -  s)  (1  -  r) 

r,  s  being  local  coordinates  as  indicated  in  Fig.  4.2. 

There  is  a  one-to-one  correspondence  between  points  in  the 
square  in  Fig.  4.2  and  the  eight-node  element  shown  in  Fig.  4.1.  It 
is  the  interpolation  used  in  Eqs.  (4.1)  and  (4.2)  that  allows  the  use 
of  curved  elements  such  as  the  ones  shown  in  Fig.  4.1. 
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Fig.  4.2  Local  system  of  coordinates  for  the  eight-node  element. 
If  small  displacements  are  assumed,  the  strains  at  any  point 


in  the  body  are  given  by 
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.  3  u 


_  3  V 


3x’  3  y 


and 


'  xy 


3u 

3  y 


3  v 

3  x 


(4.3) 


and,  using  Eg.  (4.2), 


m  3  N . 

e  =  y  - J- 

x  /=1  3  x 


m  3  N . 

UJ  ’  y  £  rt  vj 


and 


'  xy 


m 

=  I 
j=l 


3  N 

J 

3  y 


uj 


ii 

3  x 


VJ 


(4.4) 


The  partial  derivatives  of  the  interpolation  functions  with 
respect  to  x  and  y  may  be  calculated  from 


ii 

3  r 


ii 

3  S 


ii 

3  X 

ii  __ 

3  x  3  s 


3  x 
3  r 
3  x 


3  N,  3  y 

w  . 

3  y  3  r 

3  Nj  3  y 
3  y  ~l~s 


In  matrix  form 
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a  r 

a  x 

a  y 

a 

h 

X 

a  Nj 

T~r 

a““r 

a" 

=  j 

/V 

a  x 

3Nj 
a  s 

a  x 

3  y 

a 

N. 

j 

3  y_ 

a  s 

y“s 

T 

L 

y 

J  is  known  as  the  Jacobian  Matrix  and  can  be  derived  as  a  function  of 
r  and  s  using  Eq.  (4.1). 


J  = 


m  a  N, 

E  - 

j=l  9r 


a  n  . 


m 

E  _ 

j=i  a  s 


xj 


xj 


m 

£ 

j=l 


m 

E 


a  n  . 


yj 


9N 


'J 


(4.5) 


0=1  a  s 

Thus,  given  the  coordinates  r  and  s  of  a  point  in  local  coordinates, 

a  n.  a  n. 

3  ^  and  y  can  found,  for  the  corresponding  point  in  the  body: 


a  x 

a  n. 

fH 

1 

r — i 

II 

a  r 

a  Nj 

a  n. 

3  y 

a  s 

(4.6) 


The  relationship  between  the  strains  at  any  point  in  the  body 
and  the  nodal  displacements  may  be  written,  in  matrix  form,  as: 
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e  =  B  U1  (4.7) 

/v 

where 


G 


:x’  Ey  ’ 


^  -  Cuj.  Vj.  u2»  v2’  •  '  *  ’  *V  vm^ 


B  = 


3  N] 

FT 


3  N, 

_ C 

3  x 


3  x 


3  Nj 

FT 


3  N 


2  . 


a  y 


3  N1  3  N;  3  N2  3  N2 

3  y  3  x  FT"  FT" 


3  \ 
3  y 


3  Nm  3  N 
_ m  _ m 

3  y  3  x 


After  the  strains,  the  stresses  can  be  calculated  (assuming  a  plane 


state  of  stress  with  oz  =  txz  =  \^z 


0)  as 


o  =  D  e 

/"W  ^  /V/ 


where 


F,  =  [ox»  V  Txy^ 


(4.8) 
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The  principle  of  minimum  potential  energy  yields  the  stiff¬ 
ness  matrix  of  an  element  (for  example,  see  Ref.  6): 

K1  =  f  B7  £  jB  dv  (4.9) 

~  element*^  ~ 
volume 

The  integration  in  Eq.  (4.9)  is  conveniently  carried  out 
numerically.  In  this  analysis,  Gaussian  integration  is  employed. 
Nine  integration  points  have  been  used  for  the  eight-node  elements 
and  seven  integrations  points  for  the  six-node  elements.  The  inte¬ 
gration  points  in  the  eight-node  element  in  local  coordinates  are 
shown  in  Fig.  4.3. 

The  matrix  multiplication  in  Eq.  (4.9)  may  be  organized  as 

follows: 


Calling, 


I 


-I 


Fig.  4.3  Integration  points  for  the  eight-node  element 
in  local  coordinates. 


so  that. 


~  =  ^1’  ~2*  '  *  *  *  -W 

£■  rsl.  £l . Jli7 


and  using  Eq.  (4.9): 
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K1  =  E  w 
~  e,=i 


bT  db,  bT  d  b9 

~  I  ~  /ssl  ^  ~C 


bI  db,  bI  d  b9 

r^C  *+  **  X  ~  L.  **  f^C 


bt  db,  bI  d  b, 

~m  v  <*‘-'1  ~m  *+  **c 


b|  d  b 

~1  ~  ~m 


BI  D  Bm 

rsi.  ~  ~m 


BI D  Bm 

~>m  **  ~m 


(4.10) 


q  being  the  number  of  integration  points.  The  weight  of  integration 
point  i.  is  denoted  by  w^. 

The  matrix  in  Eq.  (4.10)  is  symmetric,  if  JMs  symmetric  as 
in  Eq.  (4.8). 


4.2.2  Modification  of  the  Element  Stiffness  Matrix  for 
Support  Conditions.  The  stiffness  matrix  derived  above  relates  the 
forces  at  the  nodes  of  an  element  to  the  corresponding  displacements. 
These  displacements  are  the  degrees  of  freedom.  If  a  displacement 
is  prescribed,  it  is  not  considered  a  degree  of  freedom  and  the  cor¬ 
responding  column  and  row  in  the  stiffness  matrix  J^1  are  deleted  after 
the  column  multiplied  by  the  prescribed  value  is  subtracted  from  the 
load  vector.  Thus,  for  example,  if  the  displacement  in  the  y-di recti  on 
of  node  1,  v^,  is  prescribed  equal  to  v^,  the  element  stiffness  matrix 
will  be  modified  as  follows: 
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k  k  k  k 
*1,1  1,2  K1 ,3  k1,4  '  • 

k2,l  k2,2  k2 ,3  k2,4  *  * 

k3,l  k3,2  k3 ,3  k3,4  *  * 

k4,l  k4,2  k4,3  k4,4  *  * 

•  •  •  • 

•  •  •  • 

♦  •  •  ♦ 

•  •  •  • 

k2m,lk2m,2k2m,3k2m,4  *  ' 
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Another  modification  in  the  load  vector  and  element  stiff¬ 
ness  matrix  must  be  introduced  in  this  analysis  to  take  into  account 
the  fact  that  some  degrees  of  freedom  may  be  the  same.  For  example, 
if  the  load  is  applied  to  the  specimen  shown  in  Fig.  2.1  using  rigid 
platens,  it  is  reasonable  to  assume  equal  displacements  for  all  nodes 
located  along  the  contact  edge  of  the  specimen  in  the  direction  of 
the  applied  load.  Thus,  all  degrees  of  freedom  can  be  condensed  in  a 
single  one  since  they  are  all  the  same.  The  process  of  condensation 
may  be  described  by  the  following  example.  Suppose  that  the  degrees 
of  freedom  v^  and  are  the  same  and  it  is  desired  to  condense  them 
in  a  single  degree  of  freedom,  say,  v^.  To  accomplish  this  it  is 
necessary  to  make  the  following  modifications  in  the  element  stiff¬ 
ness  matrix  and  load  vector.  The  column  and  row  corresponding  to  the 
degree  of  freedom  v^  must  be  added  to  the  column  and  row  corresponding 
to  the  degree  of  freedom  v^.  The  column  and  row  corresponding  to  the 
degree  of  freedom  v^  is  them  deleted.  Thus, 


x1 

X1 

y'Uy1 

y  2 

x1 

X2 

k.  , 

1  *  A 

(kl,2+kl,4) 

kl ,3  ’  •  ' 
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= 
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Note  that  the  remaining  stiffness  matrix  is  still  symmetric,  if  the 
original  matrix  is  symmetric. 

Since,  sometimes,  there  are  some  degrees  of  freedom  which  are 
the  same  but,  in  addition  belong  to  different  elements,  a  trick  is 
used  to  overcome  the  difficulty.  A  one-degree-of-"freedom  node  number 
is  created  in  the  program  for  each  set  of  degrees  of  freedom  to  be 
condensed.  An  element  node  number  corresponds  to  a  structure  node 
number.  For  example,  the  node  number  1  of  element  28,  see  Fig.  4.1, 
corresponds  to  node  number  87  in  the  structure.  Thus,  in  the  program, 
the  sets  of  degrees  of  freedom  to  be  condensed  are  forced  to  correspond 
to  the  created  node  number  and,  therefore,  during  the  process  of 
assembly  and  solution  of  the  system  of  equations,  only  one  degree  of 
freedom  will  be  processed. 

4.2.3  Assembly  and  Solution  of  the  Structure  System  of 
Equations.  The  assembly  and  solution  of  the  structure  system  of 
equations  is  a  straightforward  process.  A  detailed  description  of 
this  process  is  not  of  interest  in  this  study  and  therefore  only  a 
brief  presentation  will  be  given  here.  (See  Ref.  7  for  more  information). 

The  Frontal  Solution  Method  is  used.  In  this  method  the 
assembly  and  Gauss  elimination  are  performed  at  the  same  time,  whereas 
in  the  usual  solution  method  the  stiffness  matrix  for  the  structure 
is  first  assembled  and  then  Gauss  elimination  is  carried  out. 

If,  as  elements  are  processed  one  after  the  other,  a  node 
appears  for  the  last  time,  the  degrees  of  freedom  associated  with 
the  node  may  be  eliminated  and  the  corresponding  equations  are  removed 


42 


and  saved.  After  the  last  element  has  been  processed,  back-substi¬ 
tution  yields  the  nodal  displacements.  Note  that  in  the  Frontal 
Solution  Method  it  is  element  numbering  that  is  crucial  and  not  node 
numbering.  The  method  is  known  to  be  more  efficient  than  the  more 
common  band  solvers,  if  elements  with  midside  nodes  are  used  as  in 
the  present  analysis  (see  Ref.  7). 

4.2.4  Calculation  of  Stresses  and  Strains  in  the  Body.  As 
explained  in  section  4.2.1,  if  the  nodal  displacements  of  an  element 
are  known,  the  strains  and  stresses  at  any  point  within  the  element 
can  be  calculated  using  Eqs.  (4.7)  and  (4.8)  respectively.  In  this 
study,  the  calculations  are  carried  out  at  the  integration  points. 

If  the  stresses  at  all  integration  points  in  the  body  do  not 
exceed  the  material  elastic  limits,  the  linear  analysis  is  sufficient 
to  obtain  the  behavior  of  the  concrete  model.  On  the  other  hand,  if, 
at  any  integration  point,  the  stresses  exceed  the  material  elastic 
limits, a  nonlinear  analysis  must  be  performed.  This  will  be  discussed 
in  Section  4.3. 

4.2.5  Accuracy  Test  of  the  Linear  Finite  Element  Analysis. 

The  problem  depicted  in  Fig.  4. 4  was  considered  in  order  to  check  the 
accuracy  of  the  linear  finite  element  analysis.  It  involves  a  circular 
inclusion  in  a  matrix  subjected  to  uniform  tension  in  one  direction. 

An  exact  solution  of  this  problem  in  the  case  in  which  tension  is 
applied  at  an  infinite  distance  away  from  the  inclusion  (a  =  °°)  is 
available  in  the  literature  (see  Refs.  8,  9).  Details  of  the  derivation 
of  the  exact  solution  may  be  found  in  these  references.  The 


i 
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Fig.  4.4  A  test  in  the  elastic  analysis. 

displacements  of  points  A  and  B  were  obtained  from  the  exact  solu¬ 
tion  and  then  using  the  linear  finite  element  analysis.  As  it  is 
impossible  to  apply  the  tension  at  an  infinite  distance  away  from  the 
inclusion,  several  finite  element  analyses  were  performed  for  increas¬ 
ing  values  of  a  (see  Fig.  4.5).  The  results  for  two  different  values 
of  Poisson's  ratio  for  the  inclusion  are  shown  in  Table  4.1.  Good 
agreement  between  the  results  of  the  finite  element  analysis  and  the 
exact  solution  may  be  observed  as  the  value  of  a  is  increased.  Thus, 


1 
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Table  4.1  Comparison  between  results  of  the  finite  element 
analysis  and  the  exact  solution 


Ej  =  8.4  x  106  psi  E2  =  4.7  x  106  psi  =  0.20 


Vj  =  0.06 

EXACT  SOL. 

a 

0.67 

1.67 

3.33 

ao 

uA(in) 

2.23  x  10'3 

2.38  x  10"3 

2.46  x  10-3 

2.51  x  10'3 

Vin) 

-0.18  x  10"4 

-0.94  x  10‘4 

-1.41  x  10-4 

-1.75  x  10'4 

1 

v,  =  0.15 

1  I 

EXACT  SOL. 

a 

0.67 

1.67 

3.33 

oo 

UA<in> 

— 

2.40  x  10“3 

2.47  x  10"3 

2.52  x  10'3 

Vin) 

— 

-2.73  x  10'4 

-3.16  x  10"4 

-3.45  x  10"4 

the  finite  element  meshes  shown  in  Fig.  4.5  yield  reasonably  accurate 
resul ts. 

4.3  Nonlinear  Analysis 

The  nonlinear  analysis  is  better  understood  by  considering 
a  one-degree-of- freedom  structure.  Let  P  be  the  load  applied  to  the 
structure.  Also  let  u  be  the  degree  o*  freedom  associated  with  P. 

The  inelastic  behavior  of  this  structure  is  shown  in  Fig.  4.6.  Sup¬ 
pose  that  for  some  value  of  the  load  P  the  corresponding  displacement 
u  is  known.  It  is  desired  to  find  the  increment  of  displacement  £u 
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Fig.  4.6  P  vs.  u  diagram  for  a  one-degree-of-freedom  structure 

associated  with  the  increment  of  load  AP.  If,  for  any  given  displace¬ 
ment,  it  is  possible  to  find  the  corresponding  load  for  structural 
equilibrium  and  if  the  initial  elastic  stiffness  K  is  known,  the 
problem  may  be  solved  iteratively  as  follows: 

Step  1.  Assuming  that  the  structure  is  elastic,  the 
increment  of  displacement  Au^  associated  with  AP  is  found: 
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Step  2.  For  the  increment  of  displacement  Au^,  the 
corresponding  increment  of  load  AP^  for  equilibrium  is  found: 

Step  3.  The  residual  increment  of  load  ( AP  -  APj)  is  applied, 
assuming  again  elastic  behavior,  and  the  increment  of  displacement 
Au^  is  calculated: 

u  (AP  -  APj) 

A  2  i< 

Now  repeat  steps  2  and  3  until  convergence  is  obtained  and  the  value 
m 

of  Au  =  I  Au.  is  determined,  n  being  the  number  of  iterations. 
i=l  1 

The  procedure  described  above  is  usually  referred  to  as  the 
"Constant  Stiffness  Method".  It  is  used  in  a  similar  way  in  the  non¬ 
linear  analysis  of  the  concrete  model  shown  in  Fig.  2.1.  Suppose  that 
for  some  load£  on  the  model  the  nodal  displacements  U^  as  well  as  the 
stresses  and  strains  £  and  £at  all  integration  points  are  known.  It 
is  desired  to  find  the  increments  of  nodal  displacements  AU  and  the 
increments  of  stresses  and  strains  Aa  and  Ac  corresponding  to  some 
increment  of  load  AP.  The  problem  is  solved  as  follows: 

Step  1.  Assuming  that  the  structure  is  elastic,  a  linear 
analysis  (as  described  in  Section  4.2)  is  performed  and  the  increments 
of  nodal  displacements  AU,  associated  with  the  load  increment  AP  are 
found: 

jSi  ■  tsr1  a 
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Step  2.  For  each  element  and  at  each  integration  point: 

Using  AUj  and  Eq.  (4.7)  ACj  is  found.  The  increments  of  stresses 
using  the  elastic  material  properties  (Ao^e=  _D  Ae^,  see  Eq.  (4.8))  are 
calculated.  Using  the  constitutive  equations  developed  in  Chapter  3, 
the  increments  of  stresses  AOj  corresponding  to  Acj  are  evaluated.  The 
residual  increments  of  stresses  (Aoje  -  Ao^)  are  determined  and  then 
for  each  element,  the  corresponding  nodal  residual  forces  are  calcu¬ 
lated.  Applying  the  residual  forces,  another  linear  analysis  is 
performed  and  the  increments  of  nodal  displacements  AU2  are  "Found. 

Repeat  steps  1  and  2  until  a  convergence  criterion  is  satisfied. 

The  increments  of  nodal  displacements,  the  increments  of 
strains  and  the  increments  of  stresses  are  given  by: 

n 

AU  =  z  AU. 

~  1=1  ~'1 

n 

Ac  *  Z  Ac,- 

~  i=l  ~1 

n 

Ao  *  E  Ao- 

~  i=l 

n  being  the  number  of  iterations. 

4.3.1  Implementation  of  the  Constitutive  Model  in  the  Non¬ 
linear  Analysis.  As  it  may  be  concluded  from  the  above  description  of 
the  nonlinear  analysis,  a  very  important  part  is  the  calculation  of 
the  increments  of  the  stresses  on  the  material  from  the  increments 
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of  the  strains.  The  constitutive  model  described  in  Chapter  3  will 
be  used  here  to  accomplish  this  task. 

It  is  evident  from  the  description  of  the  constitutive  model 
that  the  behavior  of  the  material  deoends  on  the  loading  path.  The 
state  of  the  material  is  defined  by  the  stresses  _o  (a  plane  state  of 
stress  is  considered  so  that  C33  =  133  =  ^  =  0): 


=  [Ojj,  a22,  t12] 


(4.11) 


and  the  hardening  parameter  k.  Oj.,  o22  are  the  normal  stresses  and 


x12  is  the  shear  stress.  The  hydrostatic  stress  is 


=  3  (on  +  a22) 


and  the  shearing  stress  intensity  is  given  by 


(4.12) 


T  - 
1  "  6 


/  \2  .  2  ,  2  .  c  2 
(on  -  o22)  +  on  +  °22  +  b  T 


12 


1/2 


(4.13) 


The  assumed  yield  function  F  is 


F  (o,  k)  =  f  (o,T)  -  k. 


(4.14) 


For  the  state  of  stress  F  (a,  k)  5  0.  The  implementation 
of  the  constitutive  model  in  the  program  can  be  described  as  follows: 


1 
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First,  elastic  increments  of  stress  are  calculated  from  the 
increments  of  strain  Ac  using  Eq.  (4.8): 


Ao 


6  6  6 
Aon  ,  AO£2  >  Ax  ^2 


1  T 


Second,  the  yield  function  F  is  calculated  at  the  state  of 
stress  a.  If  F  (c,  k)  <  0,  the  behavior  is  elastic  if  the  increments 
of  stress  a oe  is  not  large  enough  for  a  +  Aae  to  be  beyond  the  yield 
surface,  i.e.,  if  F  (  a  +  Aae,  k)  <  0.  Thus 


Ac  =  Ce  A oe 


(4.15) 


C_e  being  the  elastic  compliance  matrix  (symmetric,  3x3)  and 

T 

=  [£H’  e22’  yl 2  ^ 

T 

~e  =  ^£11’  a£22*  hyl2^ 

(eu,  e22  are  norma^  strains  and  y^  is  the  shear  strain.)  The 
matrix  C?  is  given  by 


-v 


-v  0 
1  0 


0  0  2(l+v) 


(E  and  v  are  Young's  modulus  and  Poisson's  ratio  respectively.) 

The  normal  strain  and  the  increment  of  normal  strain  in  the  out-of¬ 
plane  direction  are  given  by 
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C33  =  T^v  {ell  +  Z22] 

'C32  =  Tv  ^A£ll  +  Ae22^ 
Eq.  (4.15)  may  be  rewritten  as 


(4.16a) 

(4.16b) 


A o6  =  De  Ae  (4.17) 

G  6  - 1 

where JD  =  [£  ]  is  the  elastic  rigidity  matrix  (symmetric,  3  x  3). 
In  the  case  that  F  (a  +  Ao  ,  k)  >  0,  the  increments  of  stress 
must  be  divided  into  two  parts:  one  to  bring  the  state  of  stress  to 
the  yield  surface  (elastic  behavior  as  described  before)  and  the 
remainder  resulting  in  elastoplastic  behavior  as  described  below. 

If  F(o,  k)  =  0  (the  state  of  stress  is  on  the  yield  surface), 
the  following  cases  are  distinguished: 

1)  Loading: 

T 

Aoe>  0 


3  F 

3  o 


with 
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9  r 

3  o 


3  F 


3  F 


3  F 


3a 


11 


3a 


22 


3T 


12 


The  constitutive  equations  are 

A£  =  (C6  +  CP)  Aa  (4.18) 

where  A  a  =  [Aa^,  Ao22»  At12^ 


are  the  increments  of  stress  that  correspond  to  Ac.  cf  is  the  two- 
dimensional  plastic  compliance  matrix  (3x3  )  and  C6  is,  again,  the 
elastic  compliance  matrix.  In  order  to  calculate  the  plastic  incre¬ 
ment  of  normal  strain  in  the  out-of-plane  direction,  Ae!^,  the  neces 
sary  entries  in  the  three-dimensional  plastic  compliance  matrix  will 
also  be  computed.  Thus,  with  the  notation 


Ai 


Ai 


Bi 


Bi 


Si  1 

—  +  -  B 

2  T  3 


Si 


2  T 


3  f  Si 

3  T  2  T 

2 


3  f 

Si 

3  T 

2T 

i  =  4 

1  3  f  | 

3  3  T 

,  i  =  4 


i  =  1,  2,  3 


(4.19a) 


i  =  1,  2,  3 


(4.19b) 
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(5^  —  —  o ,  —  o,  —  o 9  ri2’  ^  and  8  are  the 

plastic  shear  modulus  and  the  inelastic  dilatancy  factor  respectively 
as  defined  in  Chapter  3.)  The  entries  of  the  three-dimensional  matrix 
cf  may  be  written  as  (see  Eq.  (3.19)): 


cPj  =  A.  Bj ,  i,  j  =  1,  2,  3,  4.  (4.20) 

This  is  a  4  x  4  matrix.  Since  Ao^  is  zero,  the  third  column  may  be 
deleted.  Keeping  rows  1,  2,  4  yields  the  two-dimensional  plastic 
compliance  matrix  in  Eq.  (4.18).  Row  3  is  used  in  order  to  calculate 
the  plastic  increment  of  strain  in  the  out-of-plane  direction  Ae^. 
Eq.  (4.18)  may  be  rewritten  as 


Ao 


(4.21) 


where  £ep  =  [C6  +  cf]"1,  is  the  elastoplastic  rigidity  matrix.  In 
this  study,  the  matrix  C,p  is  not  necessarily  symmetric  and,  therefore, 
the  rigidity  matrix  £ep  is  also,  in  general  not  symmetric. 

The  increment  of  the  hardening  parameter  k  is  obtained  as 


Ak 


T 


A  a. 


(4.22) 


2)  Unloading: 


2 

3  o 


f" 

*  .  i 


T 


Ao6  <  0. 
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The  behavior  is  elastic  (Eq.  (4.17)). 
3)  Neutral  loading: 


3  F 

3  a 


Ac  =  0. 


Again,  elastic  behavior  is  obtained  (Eq.  (4.17)). 

4.3.2  Calculation  of  Residual  Forces  from  Residual  Stresses. 
After  the  residual  stresses  are  calculated  at  each  integration  point, 
the  corresponding  nodal  residual  forces  for  each  element  may  be  found 
using  the  equation: 

AF1  =  K1  AU1. 


Since 


Ac 


f  B  D  B  dv,  A.e 
voT  ~  ~  ~~ 


=  D  Ae, 


=  B  Au  and 


AF1  =  /  BT  D  B  dv  AU1  =  /  BT  D  B  Au1  dv 
vol  vol 


AF^  =  /  B^  Ac  dv 
vol 

-i  j 


(4.23) 


where  AF  is  the  vector  of  nodal  residual  forces  for  element  i  and 
Ac  represents  the  residual  stresses  in  the  element.  The  integration 
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is  carried  out  numerically: 


H  T 

E  v«  B  4o. 
£=1  ~ 


(4.24) 


Ao^  is  the  vector  of  residual  stresses  at  integration  point  Z. 

4.3.3  Convergence  Criterion.  In  this  study  the  convergence 
test  is  based  on  a  relative  error  estimate  calculated  from  the  nodal 
displacements  and  increments  of  nodal  displacements  of  the  structure. 
The  ratio  of  the  square  root  of  the  sum  of  the  squares  of  the  incre¬ 
ments  of  the  nodal  displacements  to  the  square  root  of  the  sum  of  the 
squares  of  the  nodal  displacements  in  an  iteration  is  obtained.  If 
it  is  less  than  a  specified  tolerance  then  iterations  are  terminated. 
The  value  of  the  tolerance  is  chosen  so  that  displacements  are  deter¬ 
mined  accurately  and  residual  nodal  forces  are  also  relatively  small 

at  termination  of  the  iterations.  The  tolerance  value  used  in  this 
-8 

analysis  is  10  .  After  termination  of  the  iterations,  the  remaining 

vector  of  residual  nodal  forces  is  added  to  the  next  load  increment. 

4.3.4  Failure  Criteria.  As  discussed  in  Chapter  3,  the 
hardening  parameter  k  assumes  values  between  0  and  1.  k  =  0  indicates 
the  initiation  of  material  inelasticity  and  k  =  1  means  that  the 
ultimate  strength  of  the  material  is  reached.  In  this  analysis  failure 
may  occur  locally  and  globally.  If,  at  an  integration  point,  ultimate 
strength  is  reached  the  material  is  assumed  to  fail  locally  at  that 
point  and  the  stresses  are  considered  as  residual  ones.  This  does 

not  necessarily  mean  a  global  failure  of  the  specimen.  It  may  still 


carry  load.  However,  if  equilibrium  cannot  be  satisfied  after  redis¬ 
tribution  of  stresses,  as  evidenced  by  failure  to  achieve  convergence, 
the  specimen  is  assumed  to  collapse. 


CHAPTER  5 


MORTAR-AGGREGATE  INTERFACE  MODELING 


5. 1  Introduction 

The  interface  between  the  coarse  aggregate  and  the  mortar 
is,  in  general,  the  weakest  link  in  concrete.  This  is  most  clearly 
seen  in  uniaxial  compression.  Under  this  state  of  stress,  failure  of 
the  concrete  specimen  generally  starts  with  the  development  of  bond 
cracks  at  the  interface  which  propagate  through  the  mortar.  Uniaxial 
compression  tests  have  shown  less  bond  cracks  in  high  strength  con¬ 
crete  than  in  normal  strength  concrete  (see  Chapter  2).  However, 
tests  of  the  concrete  model  specimens  in  uniaxial  compression  exhibit 
configuration  of  cracks  parallel  to  the  direction  of  the  applied  load 
and  following  a  path  along  the  interface. 

Thus  it  may  be  seen  that  it  is  very  important  to  include  in  the 
analysis  a  model  to  reproduce  the  behavior  of  the  interface  between 
the  coarse  aggregate  and  the  mortar.  In  this  study  this  is  done  by 
means  of  the  special  interface  element  shown  in  Fig.  5.1.  The  internal 
nodes  1,  3  and  5  are  attached  to  an  element  that  represents  the  aggre¬ 
gate  while  the  external  nodes  2,  4  and  6  are  attached  to  an  element 
that  represents  the  mortar  (see  Fig.  4.1  ).  The  global  coordinates  x, 
y  are  the  same  for  each  pair  of  adjacent  nodes  and  the  element  thick¬ 
ness  t  is  set  equal  to  a  very  small  number. 
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Fig.  5.1  Interface  element  in  global  and  local  systems 
of  coordinates. 

5.2  Derivation  of  the  Stiffness  Matrix 
of  the  Interface  Element 

The  interface  element  is  also  an  isoparametric  element  and 
the  derivation  of  its  stiffness  matrix  follows  a  procedure  similar 
to  that  used  in  the  derivation  of  the  stiffness  matrices  of  the  other 
elements.  The  derivation  can  be  summarized  as  follows: 

Step  1.  The  displacements  at  any  point  along  the  edges  of  the 
element  are  expressed  in  terms  of  the  nodal  displacements: 

Q 

uin  ~  ^1  U1  +  ^2  u3  +  ^3  u5  (5.1a) 


I 


I 


! 


(5.1b) 


A/ 
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T 

Un  =  v1§  u2,  v2,  .  .  ug,  v6] 

N  =  [N1,  N2,  N3] 


with 


N 


N1  0  0  0 

0  Ni  0  0 

0  0  N.  0 

0  0  0  N. 


Step  2.  The  calculated  displacements  are  transformed  to  a  local 
system  of  coordinates  that  is  tangent  and  normal  to  the  interface  at 
the  particular  point  under  consideration.  This  is  done  using  an 
appropriated  rotation  matrix: 


'uL  1 
in 

cos 

e 

sin  6 

0 

0 

G 

u  in 

v1: 

-sin 

0 

cos  9 

0 

0 

vG 

in 

in 

= 

U  out 

0 

0 

cos  0 

sin 

9 

uG 

out 

.  Vout. 

0 

0 

-sin  0 

cos 

0 

vG 

.out. 

or, 
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(5.3) 


0  is  the  angle  between  the  x  axis  and  the  tangent  to  the  interface 
at  the  point  in  study,  sin  9  and  cos  9  can  be  calculated  at  any  point 
on  the  interface  as  follows: 


cos  9 


(3  x/3  r) 
A 


sin  9 


wi  th 


A  =  [  (3x/3r)  +  (3y/3rr] 


since 

3 

x  =  I  N.  x. 

i=l  1  1 

3 

y  =  £  N.  y.  , 

i=l  1  1 

it  is  seen  that 

3  x  3  3  Ni 

ITT  =  if1  3  r  xi 
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3  y  3  3  N- 

3  r  =  ^  3  r  yi 


and  y,j  being  the  global  coordinates  of  the  pair  of  nodes  i  (e.g., 
pair  1  is  the  pair  of  nodes,  1,  2,  etc.). 

Step  3.  With  the  displacements  in  local  coordinates,  the 
strains  and  stresses  are  found  in  the  tangent  and  normal  directions. 

The  strains  are  equal  to  the  difference  between  corresponding 
displacements  divided  by  the  element  thickness. 


or. 


c  = 


T  IT 


(5.4) 


e  and  Y  are  the  strains  in  the  normal  and  in  the  tangent  directions 
respectively.  Having  calculated  the  strains,  the  stresses  can  be 


obtained  as 
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a 

E 

0 

E 

T 

0 

G 

Y 

- 

- 

- 

or, 

o  =  D  _e.  (5.5) 

o  and  x  are  the  stresses  in  the  normal  and  tangent  directions  res¬ 
pectively  and  E  and  G  are  the  extensional  and  the  shear  moduli 
of  elasticity  respectively  at  the  point  under  consideration.  In  the 
present  analysis  large  values  for  E  and  G  are  used  in  order  to  keep 
the  corresponding  displacements  of  nodes  in  each  pair  practically  the 
same.  Values  of  E  and  G  cannot  be  arbitrarily  large,  if  ill- 
conditioned  calculations  are  to  be  avoided.  An  appropriate  set  of 
values  of  E,  G,  and  t  may  be  obtained  by  trial  and  error.  In  this 
study  E  and  G  are  initially  set  equal  to  1.0  x  10®  psi  at  all  inter¬ 
face  points  and  the  element  thickness  t  is  set  equal  to  1.0  x  10"*®  in. 

Step  5.  The  element  stiffness  matrix  is  calculated. 

With  eqs.  (5.2),  (5.3)  and  (5.4)  a  relationship  between 
strains  and  nodal  displacements  may  be  derived  as  follows: 

e  =  TRUG  =  TRNU 

/s*  /“o'*  «%/ 


or 


\ 
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£ 

r^f 


B  U 


with 


(5.6) 


B  =  T  R  N 


and  the  element  stiffness  matrix  is  given  by 


K  »  /  B'  D  B  dv 
~  element  ~  ~ 
volume 


Again  the  integration  is  carried  out  numerically  using  three 
Gaussian  integration  points.  Thus,  partitioning 


B  =  [B1,  B2,  B3] 

/V.  IV 


bt  =  [(b')T.  (b2)T,  (b3)T] 


with 


B1  =  T  R  N1 


yields 


L 
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3 

K  =  Ew,  t 
~  i=l  1 


where  the  summation  is  over  the  number  of  integration  points  and 
is  the  weight  of  integration  point  i. 

5.3  Bond  Strength  Between  the  Mortar 
and  the  Coarse  Aggregate 

The  interface  may  be  submitted  to  a  combination  of  shear  and 
tensile  normal  stresses  or  to  a  combination  of  shear  and  compressive 
stresses.  Figure  5.2  shows  the  bond  failure  envelope  in  the  two 
regions.  The  envelope  is  defined  by  the  three  parameters  o^,  c  and  0, 
which  are  respectively  the  bond  tensile  strength,  the  cohesion  and 
the  angle  of  internal  friction.  These  parameters  have  been  measured 
in  bond  tests  between  coarse  aggregate  and  normal -strength  mortar. 

The  following  range  of  results  have  been  reported  for  this  type  of 
material  (see  Refs.  10,  11  and  12): 

ot  =  between  200  psi  and  400  psi 
c  =  between  300  psi  and  600  psi 

0  =  between  32  degrees  and  39  degrees 


(B1)1  D  B1  (B1)1  D  B2  ( B 1 ) T  D  B3 


(b2)T  d  b1  (b2)Td  B2  (B2)T  D  B3 


(B3)T  d  b1  (b3)T  d  b2  (b3)T  d  b3 


(5.7) 


1 
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or 


Fig.  5.2  Bond  failure  envelope  in  the  o  -  x  space. 

No  information  is  available  in  the  literature  for  the  bond 
strength  between  coarse  aggregate  and  high-strength  mortar;  thus, 
in  this  study,  different  sets  of  parameters  will  be  used  to  examine 
the  importance  of  the  bond  strength  in  the  performance  of  the  concrete 
model . 
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5.4  Implementation  of  the  Bond  Properties 
in  the  Analysis 

As  described  in  Sec.  5.2,  the  stiffness  matrix  of  the  inter¬ 
face  element  is  derived  using  initially  large  values  for  the  exten- 
sional  modulus  E  and  the  shear  modulus  G  at  all  interface  points. 

Thus,  each  pair  of  nodes  has  the  corresponding  nodal  displacements 
practically  equal  and,  therefore,  complete  attachment  (stick  condition) 
between  the  coarse  aggregate  and  the  mortar  is  achieved. 

After  the  calculation  of  the  nodal  displacements,  the  strains 
and  stresses  in  the  tangent  and  normal  directions  at  any  point  within 
the  interface  element  can  be  computed.  If  the  combination  of  normal 
and  shear  stresses  at  any  point  on  the  interface  exceeds  the  maximum 
allowed  by  the  failure  envelope  in  Fig.  5.2,  it  is  assumed  that  the 
bond  is  damaged  at  that  point  and  its  properties  are  no  longer  the 
same.  The  three  integration  points  used  in  the  derivation  of  the 
element  stiffness  matrix  are  also  used  here  as  the  locations  to  check 
if  the  bond  has  failed.  The  procedure  can  be  summarized  as  follows: 

Step  1.  After  calculating  the  nodal  displacements,  the 
strains  and  stresses  at  each  integration  point  of  each  interface 
element  are  found  in  the  tangent  and  normal  directions  using  Eqs. 

(5.6)  and  (5.5). 

Step  2.  The  calculated  stresses  are  checked  to  see  if  they 
exceed  the  maximum  allowed  by  the  failure  envelope  in  Fig.  5.2.  Here 
two  different  modes  of  bond  failure  may  occur:  shear- tension  or  shear- 
compression  failure.  If  a  shear-tension  failure  is  reached  at  an 
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integration  point,  the  values  of  E  and  6  at  that  point  are  set  equal 
to  zero  and  the  normal  and  the  tangent  stresses  become  residual 
stresses.  If  a  shear-compression  failure  is  achieved,  only  the  value 
of  the  shear  modulus  G  is  set  equal  to  zero  at  that  integration  point 
and  only  the  tangent  stress  becomes  residual  stress. 

Step  3.  Nodal  residual  forces  from  the 
residual  stresses  (see  Section  4.3.2)  and  the  stiffness  matrices  of 
the  interface  elements  in  which  the  bond  was  damaged  are  obtained, 
now  using  the  new  values  of  G  and  E  for  each  integration  point. 

Step  4.  A  new  analysis  is  performed  applying  the  residual 
forces  and  using  the  new  stiffness  matrices  for  the  damaged  interface 
elements.  If  no  more  integration  points  on  the  surface  reach  the 
failure  envelope  the  analysis  is  complete.  Otherwise,  iterations 
are  continued  until  no  further  bond  damage  occurs. 


CHAPTER  6 


ANALYTICAL  RESULTS  AND  COMPARISON 


6.1  Introduction 

The  presentation  of  the  analytical  results  in  this  chapter 
is  divided  into  two  parts.  The  results  shown  in  the  first  part  were 
obtained  using  all  the  assumptions  and  parameters  given  in  the 
preceding  chapters  (with  the  exception  of  the  bond  properties  between 
the  mortar  and  the  aggregate  which  are  discussed  in  the  next  section). 
The  second  part  of  this  chapter  shows  some  modifications  introduced 
in  the  analysis  in  order  to  obtain  better  agreement  with  the  experi¬ 
mental  results. 

6.2  Estimation  of  the  Bond  Properties 
of  the  Interface 

As  shown  in  Section  5.3,  the  bond  properties  of  the  interface 
can  be  defined  by  three  parameters,  ot,  c  and  0,  which  are,  respect- 
vely,  the  bond  tensile  strength,  the  cohesion  and  the  angle  of 
internal  friction.  Since  there  appears  to  be  no  quantitative  infor¬ 
mation  on  the  properties  of  the  bond  between  high-strength  mortar 
and  aggregates,  the  values  used  in  this  analysis  were  estimated  by 
calibration  with  experimental  results. 

Figure  6.1  shows  a  comparison  between  analytical  and  experi¬ 
mental  results  using  several  different  sets  of  bond  parameters  in 
the  case  of  uniaxial  compression.  Note  the  slope  discontinuity  in 
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the  analytical  stress-strain  curves  corresponding  to  =  300  psi, 
c  =  500  psi,  and  a^.  =  1000  psi,  c  =  1000  psi,  due  to  failure  in  the 
interface  elements.  The  value  of  the  tensile  bond  strength,  r^,  was 
assumed  to  be  equal  to  the  tensile  strength  of  the  mortar  (about  1000 
psi)  and  the  cohesion  value  was  then  adjusted  in  order  to  get  the 
same  level  of  strength  of  the  experimental  results.  A  value  of  2000 
psi  was  then  assigned  to  the  cohesion  and  the  angle  of  internal 
friction  was  kept  equal  to  35  degrees.  These  values  were  used  in  all 
other  load  cases. 

A  map  of  bond  damage  is  shown  in  Fig.  6.2  for  the  case  of 
uniaxial  compression  with  normal-strength  concrete  bond  properties. 

A  dashed  line  means  that  a  shear-compression  failure  has  occurred 
at  the  nearby  integration  point  of  the  interface  element  and  a  con¬ 
tinuous  line  means  that  a  shear- tension  failure  has  occurred  in  that 
region.  This  stage  of  loading  corresponds  to  the  points  of  slope 
discontinuity  shown  in  the  stress-strain  curves  in  Fig.  6.1. 

6.3  Analysis  Prediction  of  the  Concrete 
Model  Behavior 

The  concrete  model  behavior,  as  predicted  by  the  finite  element 
analysis,  is  presented  in  this  section.  The  elastic  and  plastic 
material  properties  used  in  the  analysis  were  obtained  from  uniaxial 
compressive  tests  on  cylindrical  specimens  as  described  in  the  prev¬ 
ious  chapters. 

Four  different  load  cases  are  studied  using  the  same  type  of 
mortar  and  coarse  aggregate:  uniaxial  compression  (o£ /a^  =  0.0), 
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6.1  Comparison  between  analytical  and  experimental  results 
using  different  bond  properties  for  the  interface 
(o2/°i  =  0.0) 
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equal  biaxial  compression  (a^/o^  =  1.0),  and  two  other  biaxial  com¬ 
pression  cases  {o^/o^  =  0.50  and  c^/c^  =  0*20).  The  load  is  applied 
on  the  specimen  in  such  a  way  that  all  the  edge  nodes  have  the  same 
displacement  in  the  direction  of  the  applied  load  and  it  is  assumed 
that  there  is  no  edge  restraint  in  the  direction  perpendicular  to  the 
direction  of  the  applied  load. 

6.3.1  Elastic  Distribution  of  Stresses.  A  graphical  repre¬ 
sentation  of  the  principal  stresses  at  the  integration  points  for 
the  four  load  cases  is  shown  in  Figs.  6.3  through  6.10.  Each  stress 
is  represented  by  an  arrow.  The  direction  of  the  arrow  coincides  with 
the  direction  of  the  principal  stress  and  the  arrow  length  is  pro¬ 
portional  to  the  magnitude  of  the  principal  stress.  A  compressive 
principal  stress  is  represented  by  an  arrow  pointing  towards  the  cor¬ 
responding  integration  point  whereas  a  tensile  principal  stress  is 
represented  by  an  arrow  pointing  in  the  opposite  direction.  Two 
figures  (with  different  scales)  are  used  to  represent  the  two  principal 
stresses  for  each  load  case.  The  first  figure  contains  the  smaller 
(in  absolute  value)  principal  stress  and  the  second  figure  contains 
the  larger  one.  The  externally  applied  stress,  Oj,  is  equal  to 
1,000  psi.  The  absolute  value  of  the  stresses  at  the  integration 
points  without  arrows  in  Fig.  6.3  is  less  than  38  psi. 

Note  that  tensile  stresses  have  occurred  only  in  the  case  of 
uniaxial  compression.  With  the  exception  of  the  case  of  equal  biaxial 
compression  (oj/oj  *  1.0),  a  slight  concentration  of  stresses  may  be 
observed  in  the  regions  between  two  circles  of  aggregate  and  the 
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Fig.  6.3  Graphical  representation  of  the  smaller  principal  stress  at 
the  integration  points  s  0. 0) - 


T 


ITT 


T 


Seal*: 


Legend: 


I - 1 - 1 

0  400  800 


empress  Ion 
tension 


Fig.  6.5  Graphical  representation  of  the  smaller  principal  stress  at 
the  integration  points  (o ?/3.  =  0.20). 
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Fig.  6.6  Graphical  representation  of  the  larger  principal  stress  at 
the  integration  points  *  0.20). 
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Fig.  6.7  Graphical  representation  of  the  smaller  principal  stress  at 
the  integration  points  *  0.50). 
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Fig.  6.9  Graphical  representation  of  the  smaller  principal  stress 
at  the  integration  points  (o^/Oj  =  1.0). 
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Fig.  6.10  Graphical  representation  of  the  larger  principal  stress 
at  the  integration  points  (o ?/o 1  =  1.0). 
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directions  of  the  principal  stresses  are  almost  coincident  with  the 
directions  of  the  applied  loads.  The  equal  biaxial  loading  case 
shows  a  more  uniform  distribution  of  stress  in  magnitude.  In  this 
case  the  directions  of  the  principal  stresses  are  approximately 
normal  and  tangential  to  the  interface  which  means  that,  practically, 
the  interface  is  subjected  to  no  shear  stresses. 

6.3.2  Evolution  of  the  Damage  in  the  Concrete  Model  During 
Loading.  As  described  in  Chapter  3,  the  hardening  parameter  k  was 
chosen  to  vary  between  0  and  1.  These  values  correspond,  respec¬ 
tively,  to  the  initiation  of  inelasticity  and  ultimate  strength  of 
the  material.  Thus,  the  magnitude  of  k  may  be  used  as  an  indicator 
of  the  proximity  to  failure  or  the  degree  of  damage  at  a  point  in 
the  specimen.  Figures  6.11  through  6.22  show  a  representation  of  the 
magnitude  of  k  at  the  mortar  elements  integration  points  for  the 
four  load  cases  and  at  diffei^.it  load  stages.  A  circle  of  radius 
proportional  to  the  value  of  k  is  used  in  the  representation.  For 
the  values  of  k  less  than  0.125  the  circles  are  not  shown. 

With  the  exception  of  the  equal  biaxial  load  case  (c^/c^  = 
1.0),  it  may  be  observed  that  the  inelasticity  of  the  material  starts 
and  develops  with  more  intensity  in  the  regions  between  two  circles 
of  aggregate  in  the  direction  perpendicular  to  o^.  Note  that  in  the 
case  of  uniaxial  compression  no  integration  point  reaches  failure  at 
ultimate  strength  of  the  specimen.  This  agrees  with  the  fact  that 
the  bond  strength  is  the  most  important  factor  influencing  the 
ultimate  strength  of  the  specimen  in  this  load  case.  The  case  of 
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equal  biaxial  compression  shows  a  more  uniform  distribution  of 
plasticity  throughout  the  specimen. 

6.3.3  Comparison  between  Analytical  and  Experimental  Results. 
Figures  6.23  through  6.26  show  the  stress-strain  curves  for  the  four 
load  cases  obtained  from  experiments  and  the  analysis.  The  average 
strains  in  the  directions  1  and  2  were  obtained  by  dividing  the  dis¬ 
placements  of  nodes  10  and  77  respectively  (Fig.  4.1)  by  the  dis¬ 
tance  from  these  nodes  to  node  1  (1.875  in.).  Nodes  10  and  77  cor¬ 
respond  to  the  points  at  which  average  strains  were  obtained  in  the 
experiments. 

For  all  the  load  cases,  the  a'alytical  results  show  less 
stiffness  than  the  experimental  ones  even  at  low  levels  of  load 
for  which  the  specimen  is  still  elastic. 

Bond  failure  for  the  load  cases  =  0.0  and  = 

was  observed  in  the  analysis  just  prior  to  failure  of  the  specimen. 

For  the  load  cases  o^/o  =  0.50  and  =  1*0  the  interface 

remained  undamaged  up  to  failure  of  the  specimen  and  the  analytical 
results  showed  higher  ultimate  'trength  and  ultimate  strains  than  the 
experimental  ones.  The  analytical  results  in  Figs.  6.23  to  6.26  were 
obtained  using  the  plastic  shear  modulus  H  as  a  function  of  the 
hardening  parameter  k  and  the  hydrostatic  stress  a  (see  Table  3.3). 

Figures  6.27  through  6.30  show  the  analytical  results  obtained 
with  the  plastic  shear  modulus  H  as  a  function  of  the  hardening  para¬ 
meter  k  only.  With  this  change  in  the  plastic  properties  of  the 
mortar,  the  load  case  o  /o^  =  0.20  did  not  show  a  bond  failure  at  the 
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Fig.  6.28  Comparison  between  analytical  and  experimental  results 
(cyoj  *  0.20;  H  =  f(k)). 


.  00  2.  00 


103 


STRAIN  *10'" 

Fig.  6.30  Comparison  between  analytical  and  experimental  results 
(< y°j  =  1.0;  H  =  f(k}). 
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interface  and  its  analytical  ultimate  strength  was  increased.  The 
load  case  c^/o^  =  0.50  showed  a  slight  decrease  in  strength  and 
less  plastic  deformation  may  be  observed  in  all  load  cases. 

The  difference  in  stiffness  between  the  analytical  and  the 
experimental  results  is  attributed  to  the  boundary  conditions  at  the 
edges  of  the  specimen  in  the  experiments.  The  next  part  of  this 
chapter  will  present  analytical  results  for  specimens  with  edge 
restraints  imposed  in  the  direction  perpendicular  to  the  direction 
of  the  applied  load. 

6.4  Analysis  of  Constrained  Specimens 

The  question  as  to  what  the  boundary  conditions  are  at  the 
edges  of  the  specimen  is  not  easily  answered  because  no  information 
is  available  about  the  degree  and  distribution  of  restraint  provided 
by  the  loading  device.  If  one  tries,  at  what  would  be  considered  as 
extreme  boundary  conditions,  to  set  equal  to  zero  the  edge  displace¬ 
ments  in  the  direction  perpendicular  to  the  direction  of  the  applied 
load,  a  singularity  in  the  displacement  field  is  obtained  at  the 
corner.  From  one  side  the  corner  must  have  the  same  displacement  as 
the  other  edge  nodes  and  from  the  other  side  this  same  displacement 
is  set  equal  to  zero.  If  the  displacements  of  the  corner  node  are 
not  prescribed  equal  to  zero,  a  concentration  of  strains  and  therefore 
stresses  is  obtained  in  the  elements  in  the  neighborhood  of  the 
corner  and  premature  failure  of  the  specimen  is  predicted. 

Any  reasonable  distribution  of  edge  restraint  is  best  justi¬ 
fied  if  good  agreement  with  the  experimental  results  is  obtained. 
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The  concrete  model  was  analyzed  with  rectangular  elements  attached 
to  the  boundary  mortar  elements  in  an  attempt  to  simulate  the 
loading  device  (Fig.  6.31).  The  elastic  properties  of  these  ele¬ 
ments  were  set  equal  to  the  elastic  properties  of  typical  steel 
(E  =  2.90  x  107  psi,  v  =  0.30)  and  the  dimension  of  the  element  in 
the  direction  of  the  applied  load  was  set  equal  to  the  dimension  of 
the  loading  device  (0.75  in.).  No  relative  displacement  was 
allowed  between  the  mortar  and  the  steel  elements  and  the  load  was 
applied  at  the  edge  nodes  of  the  steel  elements.  The  displacements 
of  the  edge  nodes  of  the  steel  elements  in  the  direction  perpendicu¬ 
lar  to  the  direction  of  the  applied  load  were  set  equal  to  zero. 

The  stress-strain  curves  obtained  with  these  boundary  con¬ 
ditions  and  the  stress-strain  curves  obtained  with  no  edge  restraint 
of  the  concrete  model  are  compared  with  the  experimental  results 
in  Figs.  6.32  through  6.35.  In  uniaxial  compression,  the  edge  con¬ 
dition  has  negligible  influence  on  the  analytical  results.  The  two 
curves  are  almost  identical.  In  the  other  three  load  cases,  the 
experimental  results  lie  in  between  the  two  sets  of  analytical  results. 
Obviously,  a  certain  degree  of  relative  displacement  has  occurred 
between  the  specimen  and  the  loading  device  in  the  experiment. 

Another  way  to  impose  edge  restraint  on  the  specimen  was 
examined  using  the  three-node  bar  element  shown  in  Fig.  6.36.  Such 
elements  were  attached  to  the  boundary  mortar  elements  in  order  to 
decrease  the  edge  displacements  in  the  direction  perpendicular  to 
the  direction  of  the  applied  load.  The  stiffness  matrix  of  this 


1 


E Xr  ER.'MENTfi'.  RESULTS 


1 


o  o  o 
o  o  o 
o  o  c 


Fig.  6.35  Comparison  between  analytical  results  with  and  without 
edge  restraints  and  the  experimental  results  (^/a^i.o) 
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Fig.  6.36  Three-node  bar  element  used  in  the  constrained  analysis 
of  the  specimen. 

element  was  obtained  numerically  using  three  Gaussian  integration 
points.  The  element  axial  stiffness  was  assumed  to  vary  linearly 
along  the  element  length.  Thus,  during  the  integration  of  the  element 
stiffness  matrix,  each  integration  point  had  a  different  value  of  the 
stiffness  assigned.  In  order  to  avoid  the  high  concentration  of 
strains  and  stresses  in  the  corner  elements,  this  axial  stiffness  was 
distributed  linearly  from  zero  at  the  corner  node  (node  100  in  Fig. 
4.1)  to  a  maximum  value  at  the  middle  of  a  side  of  the  specimen 
(nodes  13  and  94  in  Fig.  4.1).  This  maximum  value  of  axial  stiff¬ 
ness  was  calibrated  with  the  experimental  results  in  equal  biaxial 
compression  in  order  to  have  the  same  initial  elastic  stiffness  in 
the  stress-strain  curve.  A  value  equal  to  4.2  x  1071  b/tn  was  obtained, 
then  used  for  all  the  load  cases.  The  analytical  results  obtained 
under  these  assumptions  regarding  the  boundary  conditions  are  compared 
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to  the  experimental  results  in  Figs.  6.37  through  6.40.  A  better 
agreement  in  stiffness  and  ultimate  strength  may  be  observed. 

A  representation  of  the  magnitude  of  the  hardening  parameter 
k  at  different  levels  of  load  for  the  four  load  cases  is  shown  in 
Figs.  6.41  through  6.52  for  the  above  described  boundary  conditions. 
With  the  exception  of  the  uniaxial  load  case,  the  inelasticity  starts 
and  develops  with  more  intensity  at  the  corner  of  the  specimen.  The 
failure  of  the  specimen  is  seen  to  be  premature  because  of  the  edge 
conditions. 

6.5  Analysis  Shortcomings 

The  linear  finite  element  analysis  of  the  concrete  model  is 
expected  to  predict  stiffer  behavior  than  the  exact  solution  due  to 
the  discretization.  The  error,  however,  is  not  large  for  the  dis¬ 
cretization  used  in  this  analysis  and  the  results  are  accurate  enough 
as  demonstrated  in  Section  4.2.5.  This  observation  establishes  beyond 
any  doubt  that  the  specimens  were  constrained  by  the  heads  of  the 
loading  device  in  the  experiments. 

The  accuracy  of  the  nonlinear  analysis  is  mainly  influenced 
by  the  accuracy  of  the  constitutive  model  used  for  the  mortar.  The 
model  has  proved  to  be  accurate  for  normal -strength  concrete  under 
biaxial  states  of  stress  (Ref.  1)  but  has  not  yet  been  tested  for 
high-strength  mortar  under  multiaxial  states  of  stress  since  no 
experimental  data  are  available. 
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Fig.  6.42  Degree  of  damage  in  the  specimen  at  80  percent  of  the 
analytical  ultimate  strength  with  the  edge  stiffening 
elements  attached  (o2/Cj  =  0.0). 
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Fig.  6.46  Degree  of  damage  in  the  specimen  at  the  analytical 
ultimate  strength  with  the  edge  stiffening  elements 
attached  (o,/o,  *  0.20). 
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Fig.  6.49  Degree  of  damage  in  the  specimen  at  the  analytical 
ultimate  strength  with  the  edge  stiffening  elements 
attached  (c^/Oj  =  0.50). 
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Fig.  6.50  Degree  of  damage  in  the  specimen  at  58  percent  of  the 
analytical  ultimate  strength  with  the  edge  stiffening 
elements  attached  =  1*0)* 
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Fig.  6.51  Degree  of  damage  in  the  specimen  at  74  percent  of  the 
analytical  ultimate  strength  with  the  edge  stiffening 
elements  attached  (o^/a^  «  1.0). 


Fig.  6.52  Degree  of  damage  in  the  specimen  at  the  analytical 
ultimate  strength  with  the  edge  stiffening  elements 
attached  (c^/o.  =  1.0). 


The  bond  properties  used  for  the  mortar- aggregate  interface 
modeling  may  not  be  accurate  and  experimental  studies  of  the  bond 
between  the  coarse  aggregate  and  high-strength  mortar  are  yet  to  be 
performed.  An  increase  in  the  number  of  interface  points  used  to 
check  for  bond  failure  may  also  yield  more  accurate  modeling  of  the 
interface. 

The  differences  in  stiffness  between  the  analytical  and  the 
experimental  results  in  uniaxial  compression  cannot  be  attributed 
to  the  restraints  imposed  on  the  edge  of  the  specimen  only.  Other 
factors  may  be  responsible  for  the  disagreement  such  as  difference 
between  the  elastic  properties  of  the  component  materials  as  obtained 
from  tests  on  cylindrical  specimens  and  their  actual  values  in  the 
concrete  model  as  well  as  possible  eccentricities  in  the  loading 
device. 


CHAPTER  7 


SUMMARY  AND  CONCLUSIONS 


7. 1  Summary 

This  present  work  describes  a  finite  element  analysis  of  a 
high-strength  concrete  model  under  short-term  monotonic  biaxial 
compressive  loading.  The  concrete  model  consists  of  a  square  mortar 
plate  with  nine  coarse  aggregate  circular  inclusions  (Fig.  2.1). 

The  analysis  takes  into  account  the  nonlinear  behavior  of 
the  mortar  using  a  constitutive  model  proposed  in  Ref.  1.  The 
significance  of  the  bond  between  the  coarse  aggregate  and  the  mortar 
is  also  studied  using  an  interface  element  developed  in  this  work. 

The  analytical  results  are  then  compared  with  experimental  ones. 

7.2  Conclusions 

The  following  conclusions  may  be  drawn  from  this  work: 

1.  The  strength  of  the  bond  between  the  coarse  aggregate 
and  the  mortar  appears  to  be  the  most  significant  factor  influencing 
the  strength  of  the  concrete  model  in  uniaxial  compression. 

2.  Except  for  low  values  of  the  stress  ratio  02/0^  (e.g., 
°2'/ol  =  0*20),  the  analysis  predicts  that  the  bond  between  the 
coarse  aggregate  and  the  mortar  remains  intact  up  to  ultimate  strength 
of  the  concrete  and  it  is  failure  in  the  mortar  matrix  that  leads  to 
failure  of  the  model. 
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3.  Mild  stress  concentrations  were  observed  in  all  load 
cases.  Although  significant  stress  concentrations  certainly  appear 
around  aggregates  in  high-strength  concrete,  this  phenomenon, 
apparently,  is  not  captured  in  the  concrete  model  (circular  aggre¬ 
gates). 

4.  The  damage  in  the  concrete  starts  and  develops  most 
conspicuously  in  the  regions  between  two  aggregates  in  the  direction 
perpendicular  to  the  direction  of  application  of  the  largest  stress. 

5.  The  boundary  conditions  at  the  edges  of  the  specimen 
have  negligible  influence  on  the  strength  and  stiffness  of  the  con¬ 
crete  model  subjected  to  uniaxial  compression. 

6.  The  boundary  conditions  at  the  edges  of  the  specimen 
most  significantly  affect  the  strength  and  stiffness  of  the  concrete 
model  in  biaxial  compression. 
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